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Abstract In the Maxwell interaction model the collision rate is independent of the relative velocity of the 
colliding pair and, as a consequence, the collisional moments are bilinear combinations of velocity moments 
of the same or lower order. In general, however, the drift term of the Boltzmann equation couples moments 
of a given order to moments of a higher order, thus preventing the solvability of the moment hierarchy, unless 
approximate closures are introduced. On the other hand, there exist a number of states where the moment 
hierarchy can be recursively solved, the solution generally exposing non-Newtonian properties. The aim of 
this paper is to present an overview of results pertaining to some of those states, namely the planar Fourier 
flow (without and with a constant gravity field), the planar Couette flow, the force-driven Poiseuille flow, and 
the uniform shear flow. 
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1 Introduction 

As is well known, the classical kinetic theory of low-density gases began to be established as a mathematically 
sound statistical-physical theory with the work of James Clerk Maxwell (1831-1879) 1 15 , 16, 17 1. Apart from 
obtaining the velocity distribution function at equilibrium (first in 1860 and then, in a more rigorous way, 
in 1867), Maxwell derived in 1866 and 1867 the transfer equations characterizing the rate of change of any 
quantity (such as mass, momentum, or energy) which can be defined in terms of molecular properties. This 
paved the way to Ludwig Boltzmann (1844-1906) in the derivation of his celebrated equation (1872) for 
the rate of change of the velocity distribution itself. As a matter of fact, the Boltzmann equation is formally 
equivalent to Maxwell's infinite set of transfer equations. 

In his work entitled "On the Dynamical Theory of Gases," Maxwell departs from the hard-sphere model 
and writes 1.15 J7J 

"In the present paper I propose to consider the molecules of a gas, not as elastic spheres of definite 
radius, but as small bodies or groups of smaller molecules repelling one another with a force whose 
direction always passes very nearly through the centres of gravity of the molecules, and whose magni- 
tude is represented very nearly by some function of the distance of the centres of gravity. I have made 
this modification of the theory in consequence of the results of my experiments on the viscosity of air 
at different temperatures, and I have deduced from these experiments that the repulsion is inversely as 
the fifth power of the distance." 
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Maxwell realized that the hypothesis of a force inversely proportional to the fifth power of the distance, 
or, equivalently, of an interaction potential (p{r) oc r^"*, makes the nonequilibrium distribution function / 
enter the transfer equations in such a way that the transport coefficients can be evaluated without the need 
of knowing the detailed form of /. His results with this interaction model showed that the shear viscosity 
was proportional to the absolute temperature T, whereas in the case of hard spheres it is proportional to the 
square root of the temperature. As said in the preceding quotation. Maxwell himself carried out a series of 
experiments to measure the viscosity of air as a function of temperature. He writes II15II17I 

"I have found by experiment that the coefficient of viscosity in a given gas is independent of the 
density, and proportional to the absolute temperature, so that if ET be the viscosity, ET °^ p/p" 

Modern measurements of the viscosity of air show that it is approximately proportional to 7""^^ over a wide 
range of temperatures |l5]|6l, so that the actual power is intermediate between the values predicted by the 
hard-sphere and Maxwell models. 

In any case, the adoption of Maxwell's interaction model ^(r) oc r~'^ allows one to extract some useful 
information from the nonlinear Boltzmann equation, or its associated set of moment equations, for states far 
from equilibrium. Apart from its intrinsic interest, exact results derived for Maxwell molecules are important 
as benchmarks to assess approximate moment methods, model kinetic equations, or numerical algorithms 
(deterministic or stochastic) to solve the Boltzmann equation. Moreover, it turns out that many consequences 
of the Boltzmann equation for Maxwell molecules, when properly rescaled with respect to the collision fre- 
quency, can be (approximately) extrapolated to more general interaction potentials [,32] . 

The aim of this paper is to review a few examples of nonequilibrium states whose hierarchy of moment 
equations can be recursively solved for the Maxwell model. The structure and main properties of the moment 
equations for Maxwell molecules are recalled in Section |2l This is followed by a description of the solution 
for the planar Fourier flow (Section [3]), the planar Couette flow (Section |4]l, the force-driven Poiseuille flow 
(Section |5]l, and the uniform shear flow (Section |6l). The paper is closed with some concluding remarks in 
Section |2l 



2 Moment equations for Maxwell molecules 

2.1 The Boltzmann equation 

Let us consider a dilute gas made of particles of mass m interacting via a short-range pair potential (p{r). 
The relevant statistical-mechanical description of the gas is conveyed by the one -particle velocity distribution 
function /(r,v,f). The number density n(r,f), the flow velocity u(r,f), and the temperature T{r,t) are related 
to / through 

n{r,t) = J t/v/(r,v,f), (1) 



u(r,0 = ^-r /^/vv/(r,v,0, (2) 

n{r,t)kBT{r,t)=p{r,t) = -j dyVHy,r,t)f{r,y,t), (3) 

where kg is the Boltzmann constant and V(v,r,f) = v — u(r,f) is the so-called peculiar velocity. The fluxes 
of momentum and energy are measured by the pressure (or stress) tensor P and the heat flux vector q, respec- 
tively. Their expressions are 

P(r,f)=m I ^/vV(v,r,f)V(v,r,0/(r,v,0, (4) 
q(r,f) = ^/vy2(v,r,0V(v,r,0/(r,v,f). (5) 



The time evolution of the velocity distribution function is governed by the Boltzmann equation 112 11123 1 

d d 



^/ = -v.V/-|;Y^/)+W,/]- (6) 
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The first and second terms on the right-hand side represent the rate of change of / due to the free motion and 
to the possible action of an external force F, respectively. In general, such a force can be non-conservative 
and velocity-dependent, and in that case it must appear to the right of the operator d /d\. The last term on 
the right-hand side is the fundamental one. It gives the rate of change of / due to the interactions among the 
particles treated as successions of local and instantaneous binary collisions. The explicit form of the (bilinear) 
collision operator J[\\f,f] is 

(7) 



(8) 



7[v|/,/] = j dx, j dQwB{w,x) [/(v')/(v'i) -/(v)/(vi)] 

where w = v — Vi is the pre-collisional relative velocity, 

v'=v-(w-CT)a, v'l = vi 



are the post-collisional velocities, ct is a unit vector pointing on the apse line (i.e., the line joining the centers 
of the two particles at their closest approach), % is the scattering angle (i.e., the angle between the post- 
collisional relative velocity w' = v' — \\ and w), and dQ = 4|w • a\da is an element of solid angle about the 
direction of w'. Finally, B{w,x) is the differential cross section ll33t . This is the only quantity that depends on 
the choice of the potential (p{r): 



B{w,x) 



Hw,x) 

sinx 



db{w,x) 



dx 



where the impact parameter b{w,x) is obtained from inversion of 



X{w,b) 



2 dr— 

Jro(w,b) [1 



b/r 



{b/rY -4(p{r)/mw' 



,211/2 ' 



(9) 



(10) 



ro(w, b) being the distance at closest approach, which is given as the root of 1 — {b/r)^ — ^<p{r) jmw^. 

In the special case of inverse power-law (IPL) repulsive potentials of the form ^(r) = (pQ{a /r)^ , the scat- 
tering angle ;if depends on both b and w through the scaled dimensionless parameter j3 = [b / o){mw^ /2ll,(poYI ^ 
{(pQ, a, and ^ being constants), namely 



Xili)-^-2 



MP) 



dp' 



-1/2 



(11) 



where j3o(j3) is the root of the quantity enclosed by brackets. For this class of potentials the differential cross 
section has the scaling form 



2/? 



Pix) 

sinx 



dPix) 



dx 



(12) 



The hard-sphere potential is recovered in the limit ^ oo, in which case /3 = b/a, j3o(j3) — min(l,j3), 
PiX) = cos(;^/2), and B{w,x) = ^cr^ = const. On the other hand, in the case of the Maxwell potential 
(^ =4), the collision rate wB{w,x) is independent of the relative speed w, namely 

wB{wa) = Qmx), G = aV8<Po/'«. (13) 

Maxwell himself realized that if = 4 the integral in Eq. (lilt can be expressed in terms of the complete 
elliptic integral of the first kind K{z), namely 



xm = ^~2 



(2 + /34)i/4 \2 



2(2 + j34)i/2 



(14) 



Obviously, the associated function ^{x) for the IPL Maxwell potential becomes rather cumbersome and 
needs to be evaluated numerically. Nevertheless, it is sometimes convenient to depart from a strict adherence 
to the interaction potential (p{r) by directly modeling the scattering angle dependence of the differential cross 
section |[251l . In particular, in the variable soft-sphere (VSS) model one has II47II481 



/3(;c)=cos'*|, £g{x) 



-cos2('*-i)^. 



(15) 
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If one chooses t? = 1, then the scattering is assumed to be isotropic P8I and one is dealing with the variable 
hard-sphere (VHS) model [12]. However, this leads to a viscosity/self-diffusion ratio different from that of the 
true IPL Maxwell interaction. To remedy this, in the VSS model [i47i one takes ■& = 2.13986. In the context 
of me/a5?ic Maxwell models LIOJ . it is usual to take °= |g-a|~^ = l/sin(;t/2). 



2.2 Moment equations 



Let Xff{\) be a test velocity function. Multiplying both sides of the Boltzmann equation (|6| by V'^(v) and 
integrating over velocity one gets the balance equation 



where 



f = 1 Jvv/(v)/(v)=n(v/) 
is the local density of the quantity represented by V'^(v), 

= / d\\\i/{\)f{\) 

is the associated flux. 



(16) 



(17) 



(18) 



(19) 



is a source term due to the external force, and 



J^ = j dyxif{x)J[y\fJ] 

= \ jdy JdM yt/i2wB(w,i)[v/(v) + v/(vi)-v^(v')-V^(v'i)]/(v)/(vi) (20) 

is the source term due to collisions. The general balance equation (II6I 1, which is usually referred to as the 
weak form of the Boltzmann equation, is close to the approach followed by Maxwell in 1867 II17II95II . We 
say that f is a moment of order a if V'^(v) is a polynomial of degree a. In that case, if the external force 

( F) 

F is independent of velocity, the source term a^, is a moment of order a — 1 . The hierarchical structure 
of the moment equations is in general due to the flux and the collisional moment 7^,. The former is 
a moment of order a + 1, while the latter is a bilinear combination of moments of any order because of 
the velocity dependence of the collision rate wB{w,x)- In order to get a closed set of equations some kind 
of approximate closure needs to be applied. In the Hilbert and Chapman-Enskog (CE) methods 1 13 ,211 1771 
one focuses on the balance equations for the five conserved quantities, namely V'^(v) = {l,v,v^}, so that 
Jxff = {0,0,0}. Next, an expansion of the velocity distribution function in powers of the Knudsen number 
Kn = £mfp/£h (defined as the ratio between the mean free path t^^fp and the characteristic distance associated 
with the hydrodynamic gradients) provides the Navier-Stokes (NS) hydrodynamic equations and their sequels 
(Burnett, super-Burnett, . . . ). On a different vein, Grad proposed in 1949 [35.361 to expand the distribution 
function / in a complete set of orthogonal polynomials (essentially Hermite polynomials), the coefficients 
being the corresponding velocity moments. Next, this expansion is truncated by retaining terms up to a given 
order a, so the (orthogonal) moments of order higher than a are neglected and one finally gets a closed set 
of moment equations. In the usual 13 -moment approximation, the expansion includes the density n, the three 
components of the flow velocity u, the six elements of the pressure tensor P, and the three components of the 
heat flux q. The method can be augmented to twenty moments by including the seven third-order moments 
apart from the heat flux. Variants of Grad's moment method have been developed in the last few years 1791 . 
this special issue reflecting the current state of the art. 

As said above, the collisional moment J^, involves in general moments of every order. An important 
exception takes place in the case of Maxwell models, where wB{w,x) is independent of the relative speed w 
[cf. Eq. ( 112b with ^ = 4]. This implies that, if Xff{\) is a polynomial of degree a, then Jxj/ becomes a bilinear 
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combination of moments of order equal to or less than a Il381l52||95l . To be more precise, let us define the 
reduced orthogonal moments Ill9ll35ll52ll62ll95l 



Wken{c)=Nuc'L'^'^^\c^)Ytic) 



where 



c = 



2kBT 



(21) 



(22) 



is the peculiar velocity normalized with respect to the thermal velocity, ^ (c ) are generalized Laguerre 
polynomials lfT][37l (also known as Sonine polynomials), Y^(c) are spherical harmonics, c = c/c being the 
unit vector along the direction of c, and 



Nke = 



27r3/2 



1/2 



(23) 



are normalization constants, F(x) being the gamma function. The polynomials {i//a(c);Ci5 = {kijx)} form a 

complete set of orthonormal functions with respect to the inner product {F\G) = n-^^^Jdce'' F*(c)G(c). 
Let us denote as 

^ ^ dy\ifa{c)f{y) (24) 

the (reduced) moment of order a = 2k + l associated with the polynomial V^a(c) = VUfii^)- In the case of 
Maxwell models one has 



a'.a" 



(25) 



The dagger in the summation means the constraints a' + a" = a and 2 < a' < a — 2. Moreover, it is necessary 
that (£, form a triangle (i.e., \£' ~t\<e<l' + £") and ^' + ^" = ^1121. The explicit expression of the 
coefficients Cgia'a" ~ Caa"a' rather involved and can be found in Ref. fST]. A computer-aided algorithm 
to generate the collisional moments for Maxwell models is described in Ref. |66|. 

Equation ( 125b shows that the polynomials (c) are eigenf unctions of the linearized collision operator 
for Maxwell particles, the eigenvalues Xu being given by 



InnQ dxsinx^ix) 



1 + 5^5, 



worn - 



2k+i X r, f X\ 

COS ^ —Pe cos — 
2 V 2/ 



sin2*^+^'|P/.(sin 



f) 



(26) 



where Pe{x) are Legendre polynomials. Irrespective of the precise angular dependence of ^^{x), the following 
relationships hold 

hi = h+ifi, {2i+l)he = {£+l)h^u+i +au~i. (27) 
The eigenvalues can be expressed as linear combinations of the numerical coefficients 



X 

A2a =2n dx sin;t-^(Z) cos" - sin' 
Jo 2 



271 



2 

'^/3/5cos2"^sin^«ffl. 



(28) 



In particular, An = 2nQA2- The reduced eigenvalues X^f = Xu/Xi\ of order 2k + £ <6 are listed in Table [T] 
For the IPL model the coefficients must be evaluated numerically. On the other hand, by assuming Eq. il5i 
one simply has 

A2a = nm{a+-&,a+l), (29) 

where B{x, y) = r{x)r{y)/r{x+y) is the Euler beta function fT!f37'| . Table|2]gives the first few values of A2a 
for the IPL, VSS, and VHS models. It can be observed that the latter two models provide values quite close 
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Table 1 Reduced eigenvalues ^ke = ^kn/^i i for 2k + l <6 



k 


















1 





1 
i 


U 


U 





2 


3 
^ 


1 


1 


1 





3 


9 
4 


2 





1 


1 



2 
4 


7 

4 , 

7 35 /I4 
2 8/1, 


2 
1 



1 
3 
5 


3 

2 . 
11 5 A4 

4 2/l2 
<: 175 A4 

16 A, 


3 





3 

2 . 


2 
1 



2 
4 
6 


9 3 A4 
4 2 A, 
4 115 A4 
^ 16 A2 
27 357 A4 1 231 Aj 
4 16 At ^ 16 Ai 



Table 2 Coefficients (a = 1 , . . . 5) for the IPL, VSS (i5 = 2. 13986), and VHS (i5 = 1) models 



a 


ILP 


VSS 


VHS 


1 


0.685174 


0.517177 


0.523599 


2 


0.108109 


0.102913 


0.104720 


3 


0.0213745 


0.0219923 


0.0224399 


4 


0.00455921 


0.00487876 


0.00498666 


5 


0.00101072 


0.00110750 


0.00113333 



to the correct IPL ones, especially as a increases. A rather extensive table of the eigenvalues X^i for the IPL 
Maxwell model can be found in Ref. 

The most important eigenvalues are A02 and An, which provide the NS shear viscosity and thermal con- 
ductivity, namely 

= p5ij - 77ns (v,a,- + VjUi - ^ V • u5,;^ , = 

q^s = -K:KsVr, K-Ns = |^. (31) 



2.3 Solvable states 

Even in the case of Maxwell models the moment hierarchy ( 116b couples moments of order a to moments of 
order a + 1 through the flux term ^1^. Therefore, the moment equations cannot be solved in general unless 
an approximate closure is introduced. There exist, however, a few steady states where the hierarchy ( 1161) for 
Maxwell molecules can be recursively solved | 32 74 |. In those solutions one focuses on the bulk of the system 
(i.e., away from the boundary layers) and assumes that the velocity distribution function adopts a normal form, 
i.e., it depends on space only through the hydrodynamic fields («, u, and T) and their gradients. On the other 
hand, it is not necessary to invoke that the Knudsen number Kn = ^mfp/^/i (where, as said before, ^t^fp is the 
mean free path and is the characteristic distance associated with the hydrodynamic gradients) is small. 

The aim of the remainder of the paper is to review some of those solutions. All of them have the common 
features of a planar or channel geometry (gas enclosed between infinite parallel plates) and a one-dimensional 
spatial dependence along the direction orthogonal to the plates. Figure [T] sketches the four steady states to be 
considered. In the planar Fourier flow (without gravity) the recursive solvability of the moment equations is 
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(a) 



T=T' 



Hot- 



■z=+L 



9, 



Cold- 



■z=-L 



T=T 




(c) 



z=+L 



z=-L 



Fig. 1 Sketch of the steady states considered in this paper: (a) Fourier flow without (g = 0) or with (g ^ 0) gravity, (b) Couette 
flow, and (c) force-driven Poiseuille flow. 



tied to the fact that the (reduced) moments of order a = 2^ + £ are just polynomials in the Knudsen number 
(here associated with the thermal gradient) of degree a — 2 and parity £. This simple polynomial dependence 
is broken down when a gravity field is added but the problem is still solvable by a perturbation expansion in 
powers of the gravity strength. When the gas is sheared by moving plates (planar Couette flow) the dependence 
of the reduced moments on the Knudsen number associated with the thermal gradient is still polynomial, but 
with coefficients that depend on the (reduced) shear rate. In the case of the force-driven Poiseuille flow the 
nonequilibrium hydrodynamic profiles are induced by the presence of a longitudinal body force only. Again 
a perturbation expansion allows one to get those profiles, which exhibit interesting non-Newtonian features. 



3 Planar Fourier flow 

3.1 Without gravity 

In the planar Fourier flow the gas is enclosed between two infinite parallel plates kept at different temperatures 
[see Fig. [Ha)]. We assume that no external force is acting on the particles {g = 0) and consider the steady 
state of the gas with gradients along the direction normal to the plates (V zdz) and no flow velocity, 

u = 0. (32) 

Under these conditions, the Boltzmann equation ^ becomes 

vz^/(z,v)=/[v|/,/] (33) 

and the conservation laws of momentum and energy yield 

P„ = const, (34) 
q-, = const, (35) 

respectively. 

In 1979 Asmolov, Makashev, and Nosik |[9l proved that an exact normal solution of Eq. ( |33] | for Maxwell 
molecules exists with 

p = nksT = const, (36a) 
= const. (36b) 

oz 

The original paper is rather condensed and difficult to follow, so an independent derivation f/Ol is expounded 
here. 

Since in this problem the only hydrodynamic gradient is the thermal one, i.e., d^T, the obvious choice of 
hydrodynamic length is ^7- = (d^lnT)'^. As for the mean free path one can take i^fp = ^ThgT jmj Ai 1 . Both 
quantities are local and their ratio defines the relevant local Knudsen number of the problem, namely 



(37, 
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Note that An (z) n{z) °= 1 /T{z), so that ^JY{z)e{z) = const, where use has been made of Eqs. i36h Here it 
is assumed that the separation 2L between the plates is large enough as compared with the mean free path to 
identify a bulk region where the normal solution applies. Such a solution can be nondimensionalized in the 
form 



</>(c;e) = ^ 

n{z) 



2kBT{z) 



3/2 



/(z,v), 



(38) 



where c is defined by Eq. ( |22] |. All the spatial dependence of (/) is contained in its dependence on c and e. 
Thus 

df^dTdl 
dz dz dT' 



where 



-3/2 



^ ^^dTdc dTde 



Taking into account that dc/dT = —jT and de/dT = —jT 'e, one finally gets 



(39) 



(40) 



dz 



1 1 Ai 1 — 
yiksT J 2 



2+— -C + Err- 

dc de 



4>{c;e) 



(41) 



Consequently, the Boltzmann equation ( 1331) becomes 
d d 



|cj2+^-c + e^ U(c;e) = ^qj dci j [(/)(c'i;e)(/)(c';e) - (/)(ci;e)(/)(c;e)] 



/[C| (/)(£), (/)(£) 



(42) 



The orthogonal moments of (^(c; e) are 



Ht{£) = j iicv/-«(c)0(c;e), 



(43) 



so that 



Here, 



k=oe=o 



471 " 



{^ypiicjc). 



(44) 



(45) 



The moments {^k(} are a subclass of the moments defined by Eq. ( |24] |. From the definition of temperature 
and density, and the fact that the flow velocity vanishes, it follows that 



1, 



ho = ^01 = 0. 



(46) 



Taking into account the recurrence relations of the Laguerre and Legendre polynomials 101137] it is easy to 
prove that 



c~Yke{c) = {2k + i)xi/k(:{c)-2^lk { k + l+ \ ) \ifk-\,i{c), 



(47) 



where 



(Dke = i£+1) 



k + i+l 



(2£+l)(2^ + 3) 



1/2 



(Okl = {i+\) 



(2£+l)(2^ + 3) 



1/2 



(48) 



(49) 
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Multiplying both sides of Eq. (l42l) by ^fu{c) and integrating over c one gets 

d_ 
'de 



1 

^ —Xk(.'^ke + — ^ Cuvi've'-^ve^v't" ■ (50) 

The choice of the orthogonal moments {^ke} simplifies the structure of the collisional terms but, on the other 
hand, complicates the convective terms. Alternatively, one can use the non-orthogonal moments 

MM{e) = J dcc^''4(j){c;£). (51) 

In that case, one gets 

^(^2k + i-l-£-^^Mkj+i{e) = Jdc c2*^c^7[c| (£),(/>(£)] 

= Jkei£) (52) 

from Eq. (142b . Obviously, J^e is a bilinear combination of moments M^igi of order equal to or smaller than 
a^2k + t 

Equations ( 1501) and ( 1521) reveal the hierarchical character of the moment equations: moments of order a 
are coupled to moments of lower order but also to moments of order a + I. However, the hierarchy can be 
solved by following a recursive scheme. Of course, the moments of zeroth and first order, as well as one of 
the two moments of second order are known by definition [cf. Eq. ( I46H . The key point is that the moments of 
order a > 2 are polynomials in e of degree a — 2 and parity i: 

2k+i—2 

Jtu{e)= L ^r^e^ Mf'=Oifj + ^ = odd, (53) 

where ju]*^^' are pure numbers to be determined. Let us see that Eq. ( |53l l is consistent with Eq. ( |50l l. First note 
that all the moments on the left-hand side of Eq. ( l50l) have a parity different from that of but the parity is 
restored when multiplying by e or applying the operator e^d /de. Similarly, the condition 2k' + 1' + 2k" + i" = 
2k + £ assures that the product p (ji is a polynomial of the same degree and parity as those of 
Next, although the moments and are polynomials of degree 2k + t~ 1, the action of the 

operator £{2k + £ — 1 — e5 /(9e) transforms those polynomials into polynomials of degree 2fc + — 2. 

In order to complete the proof that ( 153b provides a solution to the moment equations ( 1501) . a recursive 

scheme must be devised to get the numerical coefficients \if^^ . First, notice that the left-hand side of Eq. (ISOl l 

contains at least the first power of e. So, if ^ = even one easily gets fi^ , provided that {1q is known for 

2k' < 2k + i — 2. Next, if ^ = odd, the coefficient is determined from the previous knowledge of 

H^'''^'^ for 2k' +£' <2k + £+l and of jif^'^ for 2k' + £' <2k + £- 2. Again, if£ = even and we know jU^*'^'' 

and for 2k' +£' <2k + £- 2, as well as ixf^"^ for 2k' + £' <2k + £+l,we can get and so on. As 

a starting point (moments of zeroth degree) we have = and = 0, the latter being a consequence 
of = 0. The recursive scheme is sketched in Fig. |2] where the arrows indicate the sequence followed in 

the determination of juj^^'. The open (closed) circles represent the coefficients associated with even (odd) j 
and £. Notice that the number of coefficients needed to determine a given moment is finite. In fact, the 
coefficients represented in Fig. |2] are the ones involved in the evaluation of for ^ + 2^" < 8. In practice, 
the number of coefficients actually needed is smaller since 



Aif)=0if7<maxl 1. (54) 



2k- 



3 
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2 4 6 8 10 12 14 

a=2k+l 



Fig. 2 Planar Fourier flow. Sketch of the sequence followed in the recursive determination of the numerical coefficients . 
All the coefficients of the same order a = 2^ + £ are represented by a common circle. The coefficients with j < a/3 (represented 
by circles lying below the dash-dot line) are equal to zero. So are the coefficients with j < I. 



To check the above property, note that Eq. (150b implies that 

{^t) T ( {k.e+1) (k+lI-1) (k-ll+l) (kl-l) {k-2.(.+ l) (k-ll-l)^ 

ik'i') (k"e") 




(55) 

2k'+e'+2k"+e"=2k+e, \e'-e"\<e<e'+e" 

where L.C. stands for "linear combination of". If ^j*^^^ = for j < ^{2k + £) then the first term on the right- 
hand side of Eq. ^ vanishes because j - 1 < ^{2k + 1 - 3) < ^{2k + £ - 1) < ^{2k + 1 + 1); the second 
term also vanishes because either / < ^{2k' +£') or ; - / < ^{2k" +£"). Similarly, if jif^^ = for j < £ 
then both terms on the right-hand side vanish because j — I < £ — I < £+1 and either f < £' or j — j' < £", 
respectively. The property /ij — for j < £ yields 

.^o£ = 0, £>l. (56) 

Since the non-orthogonal moments Mf^i are linear combinations of the orthogonal moments ^k'd' with 
2k' +£' <2k + £and£ + £' = even, the polynomial form (|53]l also holds for Mu: 

2k+e-2 

Mki{£)= £ mf>e^, mp =Oif j+^ = odd. (57) 

The explicit expressions for the moments of order smaller than or equal to five are given in Table |3] It also 
includes the moments of order seven, nine, and eleven corresponding to ^ = 1 and ^ = 3,5, and 7, respectively. 
The result ./#02 = or, equivalently, M02 = 5 imply that 

Pzz=P, (58) 

in agreement with Newton's law ( [30] l. The most noteworthy outcome is the linear relationship (e) = -^e 
or, equivalently, Mii(e) = — This means that Fourier's law (131b is exactly verified, i.e., 

qz=qT, (59) 

no matter how large the temperature gradient is. Likewise, ./#03 (e) = or M03 (e) = — |e imply that 

3 
5 



vl) = Uv\). (60) 
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Table 3 Planar Fourier flow. Orthogonal moments and non-orthogonal moments My {£) for 2^: + ^ < 5 and for (k,£) = 

(3, 1), (4, 1), and (5, 1). In the latter cases the numerical coefficients correspond to the IPL model. Those numerical values are, 
however, very similar to those corresponding to the VHS and VSS models |28| 
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4 ^ 28'^ 
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1 


17(225-122A4/A2) ^3 
18V35(3-2A4/A2) 


35 17(225-122A4/A2) ^3 
4 36(3-2A4/A2) 


1 



3 
5 


H 4(100-27A4/A2) ^3 


21 ^ 63225-86656A4/A2+29036(A4/A2)^ ^3 


V 5 7(1I-10A4/A2)(3-2A4/A2)'' 




4 ^ 84(11-10A4/A2)(3-2A4/A2) 
15^ 5(3125-2074A4/A2) 3 
4 84(11-10A4/A2) 


3 


1 


21.86e^-l-402.2e5 


-^e-726.2e3- 4.371 x lO^e^ 


4 


1 


-9.250e^ - 1.471 x lO^e' - 2.816 x 10'»e^ 


-2^e - 1.107 x lO^'e^ - 1.712 x 10^ - 1.436 x lO^'e'' 


5 


1 


2.038 X 10^ £5 + 1.592 x lO'e^ + 3.548 x 10''e'^ 


225225g 20286875 £3 4,593 x lO^-g^ 9.292 X lO'^e' 1.031 X lO'^e'^ 


1 (, 



iV 




0.00 0.02 0.04 0.06 0.08 0.10 



Fig. 3 Planar Fourier flow. Plot of the ratio Mux (e) /M^f (e) for, from bottom to top, = 2, 3, 4, and 5. 



The nonlinear character of the solution, however, appears through moments of fourth order and higher. Al- 
though those moments are polynomials, the velocity distribution function (l44l) involves all the powers in the 
Knudsen number e. 

Figure [3] shows the Knudsen number dependence of the moments Mj.i(e) {k = 2-5) relative to their NS 

values M^i{e) = Tfif^'^ e. It is quite apparent that, as £ increases, important deviations from the NS predic- 
tions exist, even for low values of k. This is a consequence of the rapid growth of the nonlinear coefficients of 
(£) and M^i (e) as k increases. This property also holds in the exact solution 146117111721 of the Bhatnagar- 
Gross-Krook (BGK) model kinetic equation 11111971 for the Fourier flow, where it is proven that it is closely 
related to the divergence of the CE expansion II72I . Whether or not this divergent character of the CE expan- 
sion also holds for the solution of the true Boltzmann equation for Maxwell molecules cannot be proven at 
this stage but it seems sensible to conjecture that the answer is affirmative. 
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A comparison between the analytical expressions of (e) j (e) for k = 2-5 and the numerical values 
obtained from the DSMC method for e < 0.08 shows an excellent agreement 1 28 29], thus validating both the 
practical applicability of the theoretical results in the bulk region and the reliability of the simulation method 
under highly nonequilibrium conditions. Although the validity of Fourier's law for large thermal gradients is 
in principle restricted to Maxwell models, it turns out that its practical applicability extends to other potentials, 
such as that of hard spheres |53 1. 

Before closing this Section, it is worthwhile addressing the spatial dependence of the dimensional mo- 
ments 



2kBT{z) 



= n{z) 



k+e/2 



Muieiz)). 



(61) 



Taking into account Eq. jSTj i and the fact that \jT{z)£{z) = \/ToSo = const (where Tq and eo are the temper- 
ature and the Knudsen number evaluated at z = 0, respectively), one can easily get 



MkfXz) = — eo 

m \ ml 



I 



— {«) 

"^2k+e-2-2j 



T{z) 



(62) 



where [£/2] is the integer part of £/2 and it is understood that 2k + i> 2. Equation ( |62] l shows that Mu{z) 
is just a polynomial in T{z) of degree ^ — 1 + [^/2]. As for the spatial dependence of temperature, Eq. ( I36bl l 
implies that 



T{z) = T^J\ 



2eo 



^mfp 



(0) 



(63) 



where ^^fp (0) °^ T^^'^ / p is the local mean free path at z = 0. It is convenient to define a scaled spatial variable 
i as 

dz!Xn{z!). (64) 



Since k\\{z) °= p/T(z), one has 



s = 



In terms of this new variable the temperature profile ( 163]) becomes linear. 




z-1 



(65) 



T{s) = To 1 



£0 



(66) 



Likewise, according to Eq. (162b . the moment Mj^fXs) is a polynomial in s of degree k — I + [i/2]. The exact 
solution of the BGK model for any interaction potential 14611711172 1 keeps the linear and polynomial forms 
of T{s) and Mu{s), respectively, as functions of 5. The only influence of the potential appears through the 
temperature dependence of the collision frequency (which plays the role of An), so that when applying Eq. 
( 164b one does not get Eqs. ( 1631) and ( 1651) , except for Maxwell molecules. The solution of the Boltzmann 
equation for Maxwell molecules described here can be extended to the case of gaseous mixtures [60J. 



3.2 With gravity 

Let us suppose that the planar Fourier flow is perturbed by a constant body force F = —mgi orthogonal to the 
plates and directed downwards (e.g., gravity). This situation is sketched in Fig. [Ha) with g j^O and has the 
same geometry as the Rayleigh-Benard problem. In fact, we can define a microscopic local Rayleigh number 

as 

dz A„y2W^' 
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where in the last step use has been made of Eq. dJ/l l. The dimensionless parameters y and e characterize the 
normal state of the system. Note that y = implies that either g = but e 7^ (Fourier flow without gravity) 
or e = but g j^O (equilibrium state with a pressure profile given by the barometric formula). 

The conventional Rayleigh number is Ra |y| (i/^rnfp)^- The situation depicted in Fig.[TJa) corresponds 
to a gas heated from above {d^T > ^ y > 0). If the gas is heated from below then d^T < and y < 0. We 
assume that either y > or y < but Ra < 1700, so that the gas at rest is stable and there is no convection (u = 
0). On the other hand, rarefied gases can present a Rayleigh-Benard instability under appropriate conditions 
Il20,78|. 

The stationary Boltzmann equation for the problem at hand reads 

Vz-^f{z,y)-g^nz,y)=J[y\f,f]. (68) 
In this case the moment hierarchy (fT6] l becomes 

^Mk,e+i{z)+g[2kMk-ij+i{z) + mu~i{z)] =Ju{z), (69) 

where the moments are defined by the first equality of Eq. ( 161b and Ju are the corresponding collisional 
moments. In particular, conservation of momentum and energy imply that d,MQ2 = —gn and d^M\\ = 0, 
respectively. In other words, one has 

'-^ = -SPiz), (70) 

where p = mn is the mass density, plus Eq. ( 135b . The presence of the terms proportional to g in Eqs. ( 168b and 
( 169 b complicates the problem significantly, thus preventing an exact solution (for general £ and y), even in the 
case of Maxwell molecules. On the other hand, if | y| ^ 1 one can treat it as a small parameter and carry out 
a perturbation expansion of the form 

/(z, V) = (z, v) + /(I) (z, v)y+/(2) (z, v)y2 + . . . , (71) 

M«(z) =m(?(z) +Mj^}{z)y+Ml^\z)r' + ■■■■ (72) 
Here the superscript (0) denotes the reference state without gravity, i.e., the one analyzed in Sec. 13. 11 When 
the expansion ( |72| ) is inserted into Eq. ( |69l ) one gets a recursive scheme allowing one to obtain, in principle, 

{M^^'+^'l from the previous knowledge of {Alj^f}. The technical details can be found in Ref. [83l and here 
only the relevant final results will be presented. First, it turns out that the temperature field r'^-'' of order y-i is a 
polynomial of degree j + 1 in the scaled variable 5 defined by Eq. ( 164b . This generalizes the linear relationship 

(166b found for r*^"' (i.e., in the absence of gravity). Analogously, M^^' for 2fc + ^ > 2 is a polynomial in s of 

degree k + j — 1 + [^/2] . In order to fully determine m|,^' it is necessary to make use of the collisional moments 
Jk'e' with 2k' + i' <2{2k + l + j-l). The main results are ||83] 

^T— = — - — r^y^ + ^(y^), (73) 



dz dz 5 kg 

45 
4 
"5 



^-^=^^f + ^if). (74) 



46 T 

1 + — y+^(y2) 



{vl) = \{v\ 



1 + -^7+^(7^) 



(75) 



(76) 



Equations (r73]l-(l76]l account for the first corrections to Eqs. (I36bb . (l58] l, ( 1591) , and ( |60] |, respectively, due to the 
presence of gravity. In addition, given Eqs. ( 170b and (174b . one has d^p = —pg + ^{T')- According to Eq. (175b 
the presence of gravity produces an enhancement of the heat flux with respect to its NS value when heating 
from above (y > 0). The opposite effect occurs when heating from below (y < 0). 

The above theoretical predictions for Maxwell molecules were seen to agree at a semi-quantitative level 
with DSMC results for hard spheres |82|. Moreover, an analysis similar to that of Ref. |83 1 can be carried out 
from the BGK kinetic model. The results though order y^ [85] strongly suggest the asymptotic character of 
the series (171b and ( 1721) . The theoretical asymptotic analysis of Ref. 1851 agrees well with a finite-difference 
numerical solution of the BGK equation {221 . 
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4 Planar Couette flow 

Apart from the Fourier flow, the steady planar Couette flow is perhaps the most basic nonequilibrium state. 
It corresponds to a fluid enclosed between two infinite parallel plates maintained in relative motion [see Fig. 
[Tfb)]. The plates can be kept at the same temperature or, more generally, at two different temperatures. In 
either case, in addition to a velocity profile u = Ux{z)x, a temperature profile T{z) appears in the system 
to produce the non-uniform heat flux compensating for the viscous heating. As a consequence, there are 
two main hydrodynamic lengths: the one associated with the thermal gradient (as in the Fourier flow), i.e., 
ir = {dz}nT)~\ plus the one associated with the shear rate, namely 4i = \/2k^Tjm{dux/ dz)~^ ■ Therefore, 
two independent (local) Knudsen numbers can be defined: the reduced thermal gradient defined by Eq. (137b 
and the reduced shear rate 

1 /'7'7^ 
(77) 

where for convenience here we use the collision frequency A02 associated with the viscosity, while in Eq. (|37] | 
use is made of the collision frequency X\ \ associated with the thermal conductivity. Like in the planar Fourier 
flow, the Boltzmann equation of the problem reduces to Eq. (|33l l. The conservation of momentum implies Eq. 
(|34] | as well as 

= const. (78) 



However, in contrast to Eq. (I35l l. now the energy conservation equation becomes 

-^+P,z^=0. (79) 
dz OZ 

In 1981, Makashev and Nosik r50"611 extended to the planar Couette flow the solution found in Ref. f9\ 
for the planar Fourier flow. More specifically, they proved that a normal solution to Eq. (133b for Maxwell 
molecules exists characterized by a constant pressure, i.e., Eq. (I36ab . and velocity and temperature profiles 
satisfying 

r(z)^^= const, (80a) 

OZ 

T{z)^T{z)^^ = const. (80b) 
OZ OZ 



These two equations are the counterparts to Eqs. ( [32] i and ( |36b| i. respectively. Since /I02 °= p/T, Eq. (I80ab 
implies that the dimensionless local shear rate defined by Eq. ( 177b is actually uniform across the system, i.e., 

a = const. (81) 

In the NS description, Eq. dSTb follows immediately from the constitutive equation ( [30l) and the exact conser- 
vation laws (134b and (178b . Here, however, Eq. ( 181b holds even when both Knudsen numbers £ and a are not 
small and Newton's law ( l30l) fails. This failure can be characterized by means of a nonlinear (dimensionless) 
shear viscosity V\*{a) and two (dimensionless) normal stress differences Zii 2(«) defined by 

P.^ = -^*{a)^^s^, (82) 
^^=4i(«), ^yy^ =A,{a). (83) 

p p 

Note that Eqs. ( I36ab . dTSl l, and ([SB are consistent with Eqs. (l82| and (|83]l, even though 77* (a) 7^ 1 and 
41.2(a) 7^0. 

As for Eq. (I80bb . it can also be justified at the NS level by the constitutive equation ( 131b and the exact 
conservation equation ( 179 b . Again, the validity of Eq. ( I80bl) goes beyond the scope of Fourier's law ( |3T1) . More 
specifically, Eq. ( |79b yields 

_L^ = ^r]*(ay = const. (84) 
An OZ 1 
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Equations dSObb and ( 184b are consistent with a heat flux component q, given by a modified Fourier's law of 
the form 

dT 

q, = -K*{a)K^^ — , (85) 

where K*{a) 7^ 1 is a nonlinear (dimensionless) thermal conductivity. Equations ([84l) and dSSl l allow us to 
identify the constant on the right-hand side of Eq. dSObl) . namely 



d dT 
T — T — 

dz dz 



3in{Xi\T) 



-c?6{a) = const, 



(86) 



where Q{a) = r\* {a) / K* {a) is a sort of nonlinear Prandtl number. Making use of Eq. dJTT l. Eq. ( [86l) can be 
rewritten as 



^/fe— iVfe] = — a^e{a 
aT \ J 5 



(87) 



To prove the consistency of the assumed hydrodynamic profiles ( I36al) and dSOl l we can proceed along 
similar lines as in the case of the planar Fourier flow |32,88|. First, we introduce the dimensionless velocity 
distribution function 



^(c;e,fl) = 



1 

«(z) 



2kBT(z) 



-,3/2 



/(2,V), 



where now it is important to notice that the definition ( 1221) includes the flow velocity u. Therefore, 

df^dTdl^du^dl 
dz dz dT dz dux^ 



where 



df / m \^ d^ 



du 



\2kBT I dc 



(88) 



(89) 



(90) 



The derivative df/dT is given again by Eq. ( |40] |. but now de/dT = ~T 
Eq. ( 187b . In summary. 



-1 fl. 



6c- 1^2 



£ ^a'^-d) on account of 



\2kBT J 



£ ( ^ d d \ 6 2r.r ^ - - " 

2 V^d'z^^^~- ) +7«^W37 + ^« — 



d_ 

de 



d 3 d 
de 2 dcx 



(/»(c;e,fl). (91) 



Consequently, the Boltzmann equation ( 1331) becomes 



2+ — -C + err- 

dc de 



-fl 0(a)r-- + -a- — 
5 ^ ' de 2 dcx 



(j>{c;e,a) 

■j-Q I dci J dQ -^{x) [</'(ci;e,a)(/>(c';e,a) - (/)(ci;e,a)(/)(c;e,a)] 
7[c| (/>(£,«),(/>(£, a)]. (92) 



(93) 



instead of the orthogonal moments ( 1241) . By definition, Mooo = 1, Afoio = A^oii = 0, and Mioo = j - According 



Of course, Eq. (|92] i reduces to Eq. (I42] i in the special case of the Fourier flow (a = 0). 

For simplicity, we consider now the following non-orthogonal moments of order a = 2k + i. 



Mkeh{e,a) = I dcc^^c^-''c'l^{c;e,a), 0<h<£, 



to Eq. ( 1921 ) the moment equations read 



- 2k + e-l-e^ 
de 



6 2.. X ^ 
—a 6(a)^r- 
5 ^ ' de 



Mkj+\,h + i^a {2kMk-\^i+2.h+i+hMki,h-\) 



I' 



= / dcc^''c^-''c'lJ[c\^M^Juh{£,a). 



(94) 
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While this hierarchy is much more involved than Eq. ( 152b . it can be easily checked to be consistent with 
solutions of the form 



2k+e-2 



Muh{e,a) = nij {a)e-', irij "^(a) =Oif j + £ = odd. 



(95) 



Therefore, the moments (e , fl) of order 2k + £>2 are again polynomials in the thermal Rnudsen number 
£ of degree 2k + £ — 2 and parity i. In particular, the moments of second order are independent of e, which 
yields Eqs. ^ and ^ with 



2mr\ 



Alia) =2 \m^°^^^ (a) - m''™'> (a)] , A2 (a) =3-2 [4"^^' (a) + 2m'^^^^ (a) 



The third-order moments are linear in e, as anticipated by Eq. ( 1851) with 

K (a) = m\ ' (a). 



(96) 
(97) 

(98) 



Besides, the shearing induces a component of the heat flux parallel to the flow but normal to the thermal 
gradient: 



dT 



g^-=<f'(a)K-Ns-^, (Pia) = -m\ '(a). 



(99) 



Since the reduced shear rate a is constant in the bulk domain, it is not difficult to get the spatial dependence 
of the dimensional moments 



Mmiz) = J ^/v|v-u(z)|2y-''[v,-M,(z)]V(z,v) 



n{z) 



Muh{£{z),a). 



Inserting Eq. ( [95] l into Eq. (llOOl i one has 



Mkehiz) = 



1 /t; \ k+e/2-l k-l + \l/2\ 



m 



2k+e-2~2j 



(a) Vnzjeiz) 



2k+f-2-2j 



T{z)y. 



(100) 



(101) 



As in the case of the Fourier flow, it is convenient to introduce the scaled spatial variable i through Eq. ( |64] |. 
In terms of this variable, Eqs. dSOl ) can be integrated to give 



Ux{s) = -as, T{s) = To 



1 



eo 



s — ■ 



y/2kBTo/m lOksTo 



2n/ \ 2 

-a d(a)s^ 



(102) 



where use has been made of Eqs. ( 1771) and ( 1861) . Thus, % and T are linear and quadratic functions of s, 
respectively. Analogously, \/T e is linear in s, namely 



)s. 



(103) 



Therefore, in view of Eq. (llOlb . it turns out that, if expressed in terms of s, the dimensional moment Muh{s) 
is a polynomial of degree 2k + £ — 2. Inverting Eq. ( 164b one gets the relationship between and z- 



z = Cfp(0) 



^/2k^Th/m 



So 



2 ^2kBTo/m 



a^e{a) 



lOkBTo/m 



(104) 



The solution to this cubic equation gives i as a function of z. Of course, Eqs. (1651) and (166b are recovered from 
Eqs. ( 1104b and ( 1102b . respectively, by setting a = 0. 
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Table 4 Planar Couette flow. Numerical values of the super-Burnett coefficients t]|, k'I, and ©2 = ~ '^2' according to the 
Boltzmann equation (for the IPL Maxwell potential) and several approximations 



Coefficient 


Boltzmann 


BGK model 


Ellipsoidal Statistical model 


13-moment 


R13-moment 




-3.311 


-3.60 


-4.20 


-2.60 


-3.311 


k| 


-7.260 


-6.48 


-7.88 


0.42 


-7.256 


92 


3.948 


2.88 


3.68 


-3.02 


3.945 



The key difference with respect to the Fourier flow case [cf. Eq. dSTl iI is that, while the coefficients fnf^'' 
were pure numbers, now the coefficients M^j'^'^^a), as well as 0(a), are nonlinear functions of the shear 
Knudsen number a (of parity equal to h). Unfortunately, the full dependence vn^^^^'^ (a) cannot be recursively 
obtained from Eq. ( 1941 ) because the term headed by a^d{a) couples m^l'^''^\a) to m^l'_f^^^''^\a). On the other 
hand, since that coupling is at least of order a^, it is possible to get recursively the numerical coefficients of 
the expansions of m^^^^ (a) in powers of a 

mf'\a) = l^mf"^a', mf'^=Oifi + h = odd. (105) 

!=0 

I — I 

Obviously, the coefficients m^Q are those of the Fourier flow and are obtained from the scheme of Fig. |2| 
Their knowledge allows one to get j , and so on. In general, the coefficients m^,- can be determined from 
the previous knowledge of the coefficients m^jff'''^ with 2k' +£' <2k + £ and 2k' +£' + f + i' <2k + £ + j + i. 
From the results to order one gets 02II88II 

149 

^*{a) = \+^*2a^ + &{a''), %* = -^, (106) 

K*{a) = l + Kia^ + ff{a\ ^| = (107) 

.(«)^l + e..^ + ^(.^), ,,^ 24854+135A./A. ^ ^^^3^ 

OjUU 

Ai{a) = ^a^ + &{a^), A2{a) = ^a^ + & [a^) , (109) 

^{a) = '^-a + e{a^). (110) 

The coefficients in Eqs. ( 11091 ) and dllOb are of Burnett order, while those in Eqs. ( ll06l) -( fT08] l are of super- 
Burnett order. Table |4]compares the latter coefficients with the predictions of the BGK kinetic model II14II241 
[3211451 . the Ellipsoidal Statistical (ES) kinetic model Il3lll54t first proposed by Holway ir7]|18ll40ll . Grad's 13- 
moment method II63II64II . and the regularized 13-moment (R13) method II79II80||8T1| . Despite their simplicity, 
the BGK and ES kinetic models provide reasonable values. The 13-moment method, however, gives wrong 
signs for K"| and 02- On the other hand, the more elaborate R13-moment method practically predicts the right 
results. 

The quantities ( |82| l, ( |83l l, ( |85l l, and ( 1991 ) correspond to moments of second and third order. As an illustra- 
tion of higher-order moments, let us introduce the quantities 

Rkie) = lim _ \ ' = lim — — 

Mo2l(fl) {c^c^} 

The functions Rk{£) for k= 1,2, and 3 are listed in Table|5] while Fig.|4]shows the ratio Rk{s)/R^^, where 
Rf^=Rk{0) is the NS value. 
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Table 5 Planar Couette flow. Scaled moments -/?t{e), as defined by Eq. (TTTJ, for k = 1,2, and 3 




1 |+22.41e2 

2 f +339.8e2 + 2.184 X lO^e'* 

3 ^+4.167 X 103e2 + 7.153 x 10'*e'* + 5.977 x lO^e*" 




1.6 



1.0 



0.00 



0.02 0.04 0.06 0.08 0.10 



e 



Fig. 4 Planar Couette flow. Plot of the ratio Rt{e)/Rf^ for, from bottom to top, ^: = 1, 2, and 3. 

Analogously to the case of the Fourier flow, a comparison between the theoretical expressions of the 
orthogonal moments related to Rki£) for k = 1-3 and the numerical values obtained from the DSMC method 
for e < 0.08 presents an excellent agreement II28II29I . Although the full dependence of the moments on the 
reduced shear rate a cannot be evaluated exactly from the Boltzmann equation, this goal has been achieved 
in the case of the BGK Ill4ll32ll45ll and ES II31II32II kinetic models, where one can also obtain the velocity 
distribution function (/)(c;e,a) itself. The results show that the CE expansion of the distribution function in 
powers of both £ and a is only asymptotic. The BGK and ES predictions for the coefficients T]*(fl), K*{a), 
d{a), zii 2(f )' ™d '^i'^) compare favorably well with DSMC results for both Maxwell molecules and hard 
spheres [59] . The BGK kinetic model has also been used to assess the influence of an external body force 
F = —mgi on the transport properties of the Couette flow |[84l . The results show that such an influence tends 
to decrease as the shear rate increases in the case of the rheological properties and zii 2(a), while it is 

especially important in the case of the heat flux coefficient <P[a). 



5 Force-driven Poiseuille flow 

One of the most well-known textbook examples in fluid dynamics is the Poiseuille flow l93\. It consists of 
the steady flow along a channel of constant cross section produced by a pressure difference at the distant ends 
of the channel. At least at the NS order, essentially the same type of flow can be generated by the action of 
a uniform longitudinal body force F = —mgi (e.g., gravity) instead of a longitudinal pressure gradient. This 
force-driven Poiseuille flow [see Fig. [Tic)] has received much attention from computational 14211431151116511911 
|92l|93 and theoretical E]lll|23[30l|39l|Ml|67]|76l[80l|^^ points of view. This interest has 

been mainly fueled by the fact that the force-driven Poiseuille flow provides a nice example illustrating the 
limitations of the NS description in the bulk domain (i.e., far away from the boundary layers). 
The Boltzmann equation for this problem becomes 




(112) 
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The apparent similarity to Eq. ( 168b is deceptive: although the velocity distribution function only depends on 
the coordinate (x) orthogonal to the plates, this axis is perpendicular to the force, in contrast to the case of Eq. 
, The conservation laws of momentum and energy imply 



P;,^ = const, -^ = -pg, —+P^,— =Q. 

ax ax ax 

The appropriate (dimensional) moments are defined similarly to the first line of Eq. jlOOl) . namely 



(113) 



Muh{.x) = / dy\y~vi{x)\^W-\v,^u,{x)ff{x,y), 0<h<l 



(114) 



Because of the symmetry properties of the problem, u^{x) is an even function of x and Mici>i,{x) is an even (odd) 
function of x if f — /j is even (odd). Seen as functions of g, ih{x) — u^{0) is an odd function, while Muh{x) is 
an even (odd) function if h is even (odd). The corresponding hierarchy of moment equations reads 



^ Mk^i+i,h{x) + [2kMk-u+2.h+\{^)+hMu^h-\{x)\ 



dx 



+ 8 [2kMk_ij+ij,+ i {x)+hMk^e_-i^h-\{x)] = Juh{x)- 



(115) 



Again, moments of order 2k + 1 are coupled to moments of order 2k-\-(. + \, and so their full dependence on x 
and g cannot be determined, even in the bulk domain and for Maxwell models. However, the problem can be 
solved by a recursive scheme [86] if the moments are expanded in powers of g around a reference equilibrium 
state parameterized by m,(0) = mq, p{0) = po, and T{0) = Tq: 



u,{x) =uo + u^^\x)g + uf^ {x)g' + ■■■, p{x) = + {x)g' + p(') {x)g' + • 



(116) 



T{x) = To + {x)g' + {x)g^ + . . . , M,a{x) = Mfl +Mlllg+Mgg^ + ■■■. 



(117) 



Due to symmetry reasons, Uz (x) = if j = even, p^^\x) = T^^\x) = if j = odd, and M^^^(x) = if 

h + j = odd. The functions u^'^ (x), p^-'^ {x), and T^-'' (x) are even, while m|^|(jc) has the same parity as £ — h. 
In fact, it turns out that those functions are just polynomials of degree 2j (or 2^ — 1 if i — h = odd) II86II . 
namely 



u9Hx) = 



(=1 
j 



p'-^hx] = 



(y,2/)^2i 



i=2 



2j 



(118) 



(119) 



(=2 



)=0 



where m^-^^ =Oi{i + £ — h = odd and mjl?' = 0. The numerical coefficients m1''^'^ , p^j^^'^ , r^''^'^ , and m['^^| are 
determined recursively by inserting jl 16l) -( fTT9]) into ( 11151) and equating the coefficients of the same powers 
in g and x in both sides. This yields a hierarchy of linear equations for the unknowns. This rather cumbersome 
scheme has been solved through order g^ in Ref. |[86]| . The results for the hydrodynamic profiles and the fluxes 
are 



uJx) = Ur) + 

2770 



T{x) = To 



Pix) = po 



l+C, 



f mg 
\kBTQ 



12770 Kb Tb 



( y 

KkBToJ 



Pxx=Po{ 1 - C 



ponig^ 



Pi 



Pzzix) = PQ 



(120) 
(121) 
(122) 
(123) 
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Table 6 Force-driven Poiseuille flow. Numerical values of the coefficients Cp, Cj, Cxx, C^z, and Q, according to the Boltzmann 
equation (for the IPL Maxwell potential) and several approximations 



Coefficient 


Boltzmann 


NS 


Burnett 


13-moment 


R 13-moment 


19-moment 


BGK model 


Cp 


1.2 





1.2 


1.2 


1.2 


1.2 


1.2 


Ct 


1.0153 








0.56 


0.9295 


1.04 


0.76 




6.2602 











3.36 




12.24 




6.4777 











3.413 




13.12 


Q 


0.4 





0.4 


0.4 


0.4 


0.4 


0.4 




_pW_^a Cp-Cr / mg_\ ^ 
60r)oKQTo 3 UfiTb/ 



(124) 



2 2 

qAx) = + ^{g^), qzix) = C.mgKo + ^(^^), (125) 

where po = p(0) = mpQ/kgTQ, T]o = ?7ns(0), and Kq = K"ns(0). The numerical values of the coefficients Cp, 
Ct, Cxx, C„, and Cq are shown in Table |6] which also includes the values predicted by the NS constitutive 
equations, the Burnett equations |[96l . Grad's 13-moment method ll65l . the R13-moment method II80II81II . a 
19-moment method [391, and the BGK kinetic model |87l|. 

All the coefficients Cp, Ct, Cxx, Qz, and Q vanish in the NS description |76l. At a qualitative level, the 
main correction to the NS results appears for the temperature profile. While, according to the NS equations, 
the temperature has a maximum at the mid plane x = 0, the Boltzmann equation shows that, because of the ex- 
tra quadratic term headed by Ct, the temperature actually presents a local minimum To atx = surrounded by 
two symmetric maxima Tmax at x = zLxj^ax, where Jt'niax = \/(>Ct TIqKoTq/ pq. The relative height of the maxima 
is {Ttn-dx ^ Tq) /Tq = ICt {rngx-faTaj^k^Tof- . Therefore, the temperature does not present a flat maximum at the 
middle of the channel but instead exhibits a bimodal shape with a local minimum surrounded by two symmet- 
ric maxima at a distance of a few mean free paths. This is illustrated by Fig.|5la) for g = 0.05/?q^^/Pq^^77o. The 
Fourier law is dramatically violated since in the slab enclosed by the two maxima the transverse component of 
the heat flux is parallel (rather than anti-parallel) to the thermal gradient. Other non-Newtonian effects are the 
fact that the hydrostatic pressure is not uniform across the system {Cp 7^ 0) and the existence of normal stress 

differences (Cxx 7^ 0, Q- 7^ 0), as illustrated by Fig. |5jb) for g = O-OSp^^^ / pQ^^rjQ. Moreover, there exists a 
heat flux component orthogonal to the thermal gradient {Cq 7^ 0). 

The coefficients Cp = | = 1.2 and Cq = ^ = 0.4 are already captured by the Burnett description [86 96|. 
On the other hand, the determination of Ct, Cxx, and C„ requires the consideration of, at least, super-Burnett 
contributions (e.g., d^T /dx^), as first pointed out in Ref. II86II . However, a complete determination of those 
three coefficients requires to retain super-super-Burnett terms (e.g., d'^T /dx^). In general, in order to get the 
fluxes through order g^-i, one needs to consider the CE expansion though order 2{j -h 1) in the gradients [86j ; 
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but this would also provide many extra terms of order higher than g that should be discarded. The 13- 
moment approximation 165] is able to predict, apart from the correct Burnett-order coefficients Cp and Cq, a 
non-zero value Ct = = 0.56. The R13-moment approximation significantly improves the value of Cj- (Ct — 
i§§ = 0.9295) and accounts for non-zero values of and C,, {C^^ =55= 3.36, Q,, = ^ ^ 3.413). Although 
both values are almost half the correct ones, they satisfactorily show that Q-^ is only slighter larger than Cxx, 
implying that Cw = — (Q^ — Cxx) is rather small. The more complicated 19-moment approximation |39 1 gives 
Ct" = II = 1 .04 for Maxwell molecules but the predictions for and C„ were not explicitly given in Ref. 
EH. 

The solution to the BGK equation for the plane Poiseuille flow has been explicitly obtained through order 
II87I . The results strongly suggest that the series expansion is only asymptotic, so that from a practical 
point of view one can focus on the first few terms. The results agree with the profiles (I120b -( ll25b . except that 
the numerical values of the coefficients Cj = ^ = 0.76 and, especially, Cxx = ^ = 12.24 and = ^ = 
13.12 differ from those derived from the Boltzmann equation for Maxwell molecules, as shown in Table 
|6l Interestingly enough, the BGK value of Ct agrees quite well with DSMC simulations of the Boltzmann 

equation for hard spheres [51]. It is worth mentioning that an exact, non-perturbative solution of the BGK 

3/2 1/2 

kinetic model exists for the particular value g = 2.5240/?q /pp r?n 1.3 J. 



6 Uniform shear flow 

The uniform shear flow is a time-dependent state generated by the application of Lees-Edwards boundary 
conditions II49I . which are a generalization of the conventional periodic boundary conditions employed in 
molecular dynamics of systems in equilibrium. Like the Couette flow, the uniform shear flow can be sketched 
by Fig. [ijb), except that now density and temperature are uniform and the velocity profile is strictly linear 
Il32||94||95l . The Boltzmann equation in this state becomes 

^/fcv)+v,^/(z,v) =7[v|/,/]. (126) 



Since there is no thermal gradient, the only hydrodynamic length is = ^JlkgT /m{dux/ dz)^"^ , so that the 
relevant Knudsen number is the reduced shear rate defined by Eq. (177b . which is again constant. On the 
other hand, in the absence of heat transport, the viscous heating term Px^dux/dz cannot be balanced by the 
divergence of the heat flux and so the temperature monotonically increases with time. In other words, the 
energy balance equation ( 1791) is replaced by 

3nkR dT „ dux „ 

We can still define the dimensionless distribution function (|88] | and the dimensionless moments ( 1931) , but now 
e = (no thermal gradient), so these quantities are nonlinear functions of the reduced shear rate a only. The 
moment hierarchy is formally similar to Eq. ( 1941 ). except that the first term on the left-hand side is replaced 
by a term coming from the time dependence of temperature, in the same way as the first term on the left-hand 
side of Eq. (179b is replaced by that of Eq. ( 1127b . More explicitly, the moment equations are 

^+0 E,Muh+a{2kMk-i/+2M+\ +hMujr-i) = -^luh{a), (128) 

where t, = ?LQ2^dlnT /dt. In contrast to Eq. ( |94l ), the hierarchy ( 11281) only couples moments of the same order 
2k + £ and of lower order, so it can be solved order by order. Note that moments of odd order (like the heat 
flux) vanish because of symmetry. 

Specifically, setting {k,£,h) = (1,0,0), (0,2,0), (0,2, 1), and (0,2,0) in Eq. ( [128]) , one gets 

_ , , 1 1 _ , , 1 a _ , , 1 l+3£(fl) 



where ^ (a) is the real root of the cubic equation 3^ (1 + ^ ) = 2a , namely 



4 , 
^(fl) = -sinh-" 



7Cosh-'(l+9a^) 
6 



(130) 



22 



Andres Santos 



1.2 




0.2- 

o.qI ■ ' ■ ' • ' • 

0.0 0.5 1.0 1.5 2.0 



a 

Fig. 6 Uniform shear flow. Plot of the nonlinear shear viscosity T)*(a) (solid line) and the normal stress difference Ai (a) (dashed 
line). 



30 

20 - 
10 - 




Ql . 1 , 1 , ^ , 1 

2 4 6 8 

a 

Fig. 7 Uniform shear flow. Plot of the exponent T(a) defined by Eq. (|T32) for the IPL model (solid line) and the VHS model 
(dashed hne). The dotted horizontal hne t = 2 intercepts the curves at Oc — 6.846 and ac ~ 7.746, respectively. 

Upon deriving Eqs. ( |129l l and dOOI l use has been made of the properties Mioo = -^loo = 0, /o20 = -^02 (A^020 
■^021 = -■^02-^^021. and 7022 = -^02 {M022 - j)- 

The nonlinear shear viscosity T]*(a), defined by Eq. ( [82l) . and the normal stress differences zli 2(a), defined 
by Eq. ( 183b , are now explicitly given by 
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The dependence of Tl*{a) and Ai (a) on the reduced shear rate a is shown in Fig.|6] For small shear rates, one 
has ^ (a) = fa^ + ff{a^), so that T]*(a) = 1 - ^a^ + &{ a^) [a) = 2a^ + 6'{a^). The latter two quantities 
differ from the ones in the Couette flow, Eqs. ( 11061) and ( 11091 ). 

Once the second-order moments are fully determined, one can proceed to the fourth-order moments. 
Although not related to transport properties, they provide useful information about the population of particles 
with velocities much larger than the thermal velocity. From Eq. (1128b one gets a closed set of nine linear 
equations that can be algebraically solved Il32ll73ll75l . Interestingly, these moments are well-defined only 
for shear rates smaller than a certain critical value a,, (a,, 6.846 and 7.746 for the IPL and VHS models, 
respectively). Beyond that critical value, the scaled fourth-order moments [e.g., (c^) = (V^) /{IkgT /in)^] 
monotonically increase in time without upper bound and diverge in the long-time limit. This clearly indicates 
that the reduced velocity distribution function exhibits an algebraic high- velocity tail 02ll56ll57l 

^(c;a)^c~^-'^''\ (132) 

so that those moments of order equal to or larger than 2 -|- z{a) diverge. In particular, the critical shear rate 
flc is the solution to z{a) = 2. The dependence of the exponent z{a) on the shear rate is shown in Fig.|7l The 
scenario described by Eq. (1132b has been confirmed by DSMC results for Maxwell molecules |58|. On the 
other hand, hard spheres do not present an algebraic high- velocity tail 11561 and thus all the moments are finite. 

7 Conclusions 

In this paper I have presented an overview of a few cases in which the moment equations stemming from the 
Boltzmann equation for Maxwell molecules can be used to extract exact information about strongly nonequi- 
librium properties in the bulk region of the system. This refers to situations where the Knudsen number defined 
as the ratio between the mean free path and the characteristic distance associated with the hydrodynamic gra- 
dients is finite, whilst the ratio between the mean free path and the size of the system is vanishingly small. 
In other words, the system size comprises many mean free paths but the hydrodynamic quantities can vary 
appreciably over a distance on the order of the mean free path. 

In the conventional Fourier flow problem, the relevant Knudsen number (in the above sense) is given by 
Eq. ( 137b . Regardless of the value of e, it turns out that the infinite moment hierarchy admits a solution where 
the pressure is uniform, the temperature is linear in the scaled variable i defined by Eq. ( 1641) [cf. Eq. ( I66b 1. 
and the dimensionless moments of order a are polynomials in £ of degree a — 2 [cf. Eqs. (153b and ( I57b 1. The 
latter implies that the dimensional moments of order a are polynomials in i of degree [a/2] — 1 [cf. Eq. ( I62b 1. 
In this solution the heat flux is exactly given by Fourier's law, Eq. ( |59b . even at far from equilibrium states. 
On the other hand, the situation changes when a gravity field normal to the plates is added. Apart from the 
Knudsen number e, a relevant parameter is the microscopic Rayleigh number defined by Eq. (167b . Assuming 
\Y\ ^ 1, a perturbation expansion can be carried out about the conventional Fourier flow state, the results 
exhibiting non-Newtonian effects [cf. Eq. ( 1741) 1 and a breakdown of Fourier's law [cf. d75]l1. 

When the parallel plates enclosing the gas are in relative motion, one is dealing with the so-called planar 
Couette flow. The plates can be kept at the same temperature or at different ones. In the former case the state 
reduces to that of equilibrium when the plates are at rest, while it reduces to the planar Fourier flow in the latter 
case. In either situation, the shearing produces a non-uniform heat flux to compensate (in the steady state) for 
the viscous heating, so that a temperature profile is present. This implies that two independent Knudsen 
numbers can be identified: the one associated with the flow velocity field [cf. Eq. ( |77b 1 and again the one 
related to the temperature field [cf. Eq. ( I37b 1. Now the exact solution to the moment equations is characterized 
by a constant pressure, a constant reduced shear rate a, a quadratic dependence of the temperature profile on 
the scaled spatial variable i [cf. Eq. ( I102b 1. and again a polynomial dependence of the dimensionless moments 
on the reduced thermal gradient £ [cf. Eq. (195 b 1. The dimensional moments of order a are polynomials of 
degree a — 2 in 5 [cf. Eq. dlOlM . As a consequence, Newton's law is modified by a nonlinear shear viscosity 
[cf. Eq. ( |82b 1 and two viscometric functions [cf. Eq. (1831)1. and a generalized Fourier's law holds whereby the 
heat flux is proportional to the thermal gradient but with a nonlinear thermal conductivity [cf. Eq. ( I85b 1 and a 
new coefficient related to the streamwise component [cf. Eq. (|99b1. The main difficulty lies in the fact that the 
coefficients of those polynomials are no longer pure numbers but nonlinear function of the reduced shear rate 
a. Therefore, although the structural form of the solution applies for arbitrary a, in order to get explicit results 
one needs to perform a partial CE expansion in powers of a and get the coefficients recursively. In particular, 
the super-Burnett coefficients t]| and k"| are obtained [cf. Eqs. ( 11061) and ( 11071) 1. 
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One of the nonequilibrium states more extensively studied in the last few years is the force-driven Poiseuille 
flow. A series expansion about the equilibrium state in powers of the external force allows one to get the hy- 
drodynamic profiles and the moments order by order. Calculations to second order suffice to unveil dramatic 
deviations from the NS predictions. More explicitly, the temperature profile presents a bimodal shape with a 
dimple at the center of the channel [cf. Fig. [5ja)]. Besides, strong normal-stress differences are present [cf. 



In the specific states summarized above the moment equations have a hierarchical structure since the 
divergence of the flux term, represented by V • 0^ in Eq. ( 116b , involves moments of an order higher than that 



of f . This is what happens in Eqs. (150b , ( 152b , i69i . i94i . and dl 15b . As a consequence, a recursive procedure 
is needed to obtain explicit results step by step without imposing a truncation closure. On the other hand, 
the methodology is different (and simpler!) in time-dependent states that become spatially uniform in the co- 
moving or Lagrangian frame of reference. In that case the divergence of the flux term <Pxf, involves moments 
of the same order as iff and so the hierarchical character of the moment equations is truncated in a natural way. 
This class of states include the uniform shear flow |[2l l32ll4ll l55l)57!l73l,l75P94l'951 considered in Sec. |6l and 
the uniform longitudinal flow (or homoenergetic extension) 1^,34,44,68 69 95|. In both cases the moments 
of second order (pressure tensor) obey a closed set of autonomous equations that can be algebraically solved 
to get the rheological properties as explicit nonlinear functions of the reduced velocity gradient, the solution 
exhibiting interesting non-Newtonian effects. Once the moments of second order are determined, one can 
proceed to the closed set of equations for the moments of fourth order. Here the interesting result is that those 
moments diverge beyond a certain critical value of the reduced velocity gradient |2 32 68 75], what indicates 
the existence of an algebraic high- velocity tail of the velocity distribution function li2] |32ll57l . 

Needless to say, exact results in kinetic theory are important by themselves and also as an assessment of 
simulation techniques, approximate methods, and kinetic models. In this sense, it is hoped that the states and 
results reviewed in this paper can become useful to other researchers. 
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